Turbine Collision Risk: Generate Distribution Data

Cédric Scherer https://cedricscherer.com (Leibniz Institute for Zoo and Wildlife Research)https://www.izw-berlin.de/en/home.html , Christian Voigt (Leibniz Institute for Zoo and Wildlife Research)https://www.izw-berlin.de/en/home.html
2022-03-31

Preparation

Show code
library(tidyverse)
library(here)
library(broom)
library(patchwork)

theme_set(theme_light(base_size = 11.5, base_family = "Open Sans"))
theme_update(
  panel.grid.major = element_line(size = .3, color = "grey93"),
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  strip.background = element_rect(fill = "grey60", colour = "grey60"),
  strip.text.x = element_text(hjust = 0),
  strip.text = element_text(size = 14, face = "bold"),
  axis.title.x = element_text(size = 14, margin = margin(t = 12)),
  axis.title.y = element_text(size = 14, margin = margin(r = 12)),
  axis.text = element_text(size = 12),
  legend.title = element_text(size = 14, face = "bold"),
  legend.text = element_text(size = 12),
  plot.title = element_text(size = 24, face = "bold", margin = margin(5, 0, 5, 0)),
  plot.subtitle = element_text(margin = margin(5, 0, 25, 0), size = 17),
  plot.title.position = "plot",
  plot.margin = margin(rep(12, 4))
)

## source functions to generate and visualize distributions
source(here("R", "src", "generate-distributions.R"))

Research Question

Is the predictive power of the current acoustic monitoring sufficient enough to estimate the true number of killed bats?

Predictions

The predictive power of the acoustic monitoring for the extrapolation of the expected number of stroke victims is decreasing

Open Questions

Inputs

We have the following parameters we can vary:

Simulated Data

Show code
scenarios <- 
  read_rds(here("output", "data-proc", "simulation-runs-norm.rds")) %>% 
  ## remove extremely skewed distributions
  filter(!str_detect(distribution, "_5$")) %>% 
  ## remove outdated 23% scenario
  filter(prop_monitored != .23) %>% 
  group_by(n, prop_monitored, distribution) %>% 
  mutate(diff = prop_monitored - median(prop_n_monitored)) %>% 
  ungroup() %>% 
  mutate(
    distribution = factor(
      distribution, 
      levels = c(
        "uniform", "random", 
        "inner_1", "inner_3",
        "outer_1", "outer_3", 
        "bottom_1", "bottom_3", 
        "top_1", "top_3"
      ), 
      labels = c(
        "a) uniform", "b) random",
        "c) inner — weak", "d) inner — strong",
        "e) outer — weak", "f) outer — strong",
        "g) bottom — weak", "h) bottom — strong",
        "i) top — weak", "j) top — strong"
      )
    ),
    prop_monitored_lab = paste0(prop_monitored * 100, "%"),
    prop_monitored_lab = fct_reorder(prop_monitored_lab, prop_monitored)
  )

Visualizations Bat Passes

Error Bars Deviation Bat Passes Recorded vs. Expected

Distribution (col) x Monitored Area (color) x Passes (x)

Show code
scenarios %>% 
  group_by(distribution, prop_monitored, n) %>% 
  summarize(mean = mean(prop_n_monitored), sd = sd(prop_n_monitored)) %>% 
  ungroup() %>% 
  mutate(
    base = mean - prop_monitored,
    min = base - sd, 
    max = base + sd
  ) %>% 
  ggplot(aes(factor(n), base, color = prop_monitored, group = prop_monitored)) +
  geom_hline(aes(yintercept = 0), color = "grey75", size = 1.2)  + 
  geom_line(
    aes(color = prop_monitored, 
        color = after_scale(colorspace::desaturate(colorspace::lighten(color, .6), .4))),
    size = .7, show.legend = FALSE
  ) +
  geom_pointrange(aes(ymin = min, ymax = max), size = .5) +
  facet_wrap(~distribution, nrow = 2, dir = "v") +
  scale_x_discrete(expand = c(.05, .05)) +
  scale_y_continuous(expand = c(.012, .012), breaks = -5:5 / 10) +
  # scico::scale_color_scico(
  #   palette = "bamako", end = .85, direction = -1, name = "Proportion\ncovered\nby AUD:",
  #   breaks = c(.04, seq(.05, .5, by = .05)), labels = scales::percent_format(accuracy = 1)
  # ) +
  # scale_color_viridis_c(
  #   option = "turbo", name = "Proportion\ncovered\nby AUD:",
  #   breaks = c(.04, seq(.05, .5, by = .05)), labels = scales::percent_format(accuracy = 1)
  # ) +
  scico::scale_color_scico(
    palette = "hawai", end = .85, direction = -1, name = "Proportion\ncovered\nby AUD:",
    breaks = c(.04, seq(.05, .5, by = .05)), labels = scales::percent_format(accuracy = 1)
  ) +
  guides(color = guide_legend(keywidth = unit(.6, "lines"), keyheight = unit(1.2, "lines"))) +
  labs(x = "Number of bat passes (categorical axes)", y = "Deviation from expected proportion") +
  theme(panel.spacing.x = unit(.8, "lines"), axis.text = element_text(size = 10),
        legend.text = element_text(hjust = 1))
Show code
ggsave(here("plots", "fig3_passes_recorded_mean_sd_lines_cat_alt2.pdf"), 
       width = 13.5, height = 7.3, device = cairo_pdf)
ggsave(here("plots", "fig3_passes_recorded_mean_sd_lines_cat_alt2.png"), 
       width = 13.5, height = 7.3, dpi = 800)

## continuous x axis
# scenarios %>% 
#   group_by(distribution, prop_monitored, n) %>% 
#   summarize(mean = mean(prop_n_monitored), sd = sd(prop_n_monitored)) %>% 
#   ungroup() %>% 
#   mutate(
#     base = mean - prop_monitored,
#     min = base - sd, 
#     max = base + sd
#   ) %>% 
#   ggplot(aes(n, base, color = prop_monitored, group = prop_monitored)) +
#   geom_hline(aes(yintercept = 0), color = "grey75", size = 1.2)  + 
#   geom_line(
#     aes(color = prop_monitored, 
#         color = after_scale(colorspace::desaturate(colorspace::lighten(color, .6), .4))),
#     size = .7, show.legend = FALSE
#   ) +
#   geom_pointrange(aes(ymin = min, ymax = max), size = .4) +
#   facet_wrap(~distribution, nrow = 2, dir = "v") +
#   scale_x_continuous(expand = c(.04, .04), breaks = unique(scenarios$n),
#                      guide = guide_axis(n.dodge = 2)) +
#   scale_y_continuous(expand = c(.012, .012), breaks = -5:5 / 10) +
#   scico::scale_color_scico(
#     palette = "bamako", end = .8, direction = -1, name = "Proportion\ncovered\nby AUD:",
#     breaks = seq(.05, .5, by = .05), labels = scales::percent_format(accuracy = 1)
#   ) +
#   guides(color = guide_legend(keywidth = unit(.6, "lines"), keyheight = unit(1.2, "lines"))) +
#   labs(x = "Number of bat passes", y = "Deviation from expected proportion") +
#   theme(panel.spacing.x = unit(.6, "lines"), 
#         axis.text.x = element_text(size = 9),
#         axis.text.y = element_text(size = 12),
#         legend.text = element_text(hjust = 1))
# 
# ggsave(here("plots", "figX_passes_recorded_mean_sd_lines_cont.pdf"), 
#        width = 17, height = 8.5, device = cairo_pdf)
# ggsave(here("plots", "figX_passes_recorded_mean_sd_lines_cont.png"), 
#        width = 17, height = 8.5, dpi = 800)

Boxplots Relative Difference Predicted vs. Expected

Dedicated Scenarios

Show code
colors <- scico::scico(n = 9, palette = "roma")[c(8,2)]

data_rel_scen <- 
  scenarios %>%
  filter((n == 400 & prop_monitored == .04) | (n == 100 & prop_monitored == .5)) %>% 
  mutate(
    diff_passes = ((n_monitored / prop_monitored) - n) / n,
    scenario = if_else(prop_monitored == .04, "Low", "High")
  ) %>% 
  bind_rows(
    tibble(
      distribution = factor("a) uniform", levels = levels(scenarios$distribution)),
      label = c("Overestimation", "Underestimation"),
      scenario = "High",
      diff_passes = c(1.25, -.75),
      color = colors
    )
  )

data_rel_scen %>% 
  filter(is.na(label)) %>% 
  ggplot(aes(scenario, diff_passes)) +
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf, fill = colors[1], alpha = .12) + 
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0, fill = colors[2], alpha = .12) + 
  geom_text(data = data_rel_scen, aes(label = label, color = color), hjust = .15, fontface = "bold", size = 5) +
  geom_hline(yintercept = 0, size = .6, linetype = "31", color = "grey75") +
  geom_boxplot(color = "grey45", size = .6, width = .5, 
               position = position_nudge(x = .05), 
               outlier.size = .9, outlier.shape = 1) +
  geom_point(color = "grey90", shape = "-", size = 10, 
             position = position_nudge(x = -.31)) +
  geom_point(color = "grey45", shape = "-", size = 10, alpha = .05, 
             position = position_nudge(x = -.31)) +
  stat_summary(geom = "point", shape = 18, size = 4.5, color = "#212121", 
               position = position_nudge(x = .05)) +
  facet_wrap(~distribution, nrow = 2, dir = "v") + 
  scale_y_continuous(expand = c(.04, .04), breaks = seq(-1, 3.5, by = .5)) +
  scale_color_identity() +
  labs(x = "Coverage scenario",
       y = "Relative difference of predicted versus expected bat passes") +
  theme(axis.text.x = element_text(size = 13))
Show code
ggsave(here("plots", "fig4_passes_difference_boxplots_rel_scenarios.pdf"), 
       width = 12, height = 7.7, device = cairo_pdf)
ggsave(here("plots", "fig4_passes_difference_boxplots_rel_scenarios.png"), 
       width = 12, height = 7.7, dpi = 800)

All Scenarios with 400 Bat Passes

Show code
data_rel <- 
  scenarios %>%
  filter(n == 400) %>% 
  filter(!prop_monitored == .04) %>% ## remove dedicated scenario
  mutate(diff_passes = ((n_monitored / prop_monitored) - n) / n) %>% 
  bind_rows(
    tibble(
      distribution = factor("a) uniform", levels = levels(scenarios$distribution)),
      label = c("Overestimation", "Underestimation"),
      prop_monitored_lab = factor("10%", levels = levels(scenarios$prop_monitored_lab)),
      diff_passes = c(1.25, -.67),
      color = colors
    )
  )

data_rel %>% 
  filter(is.na(label)) %>% 
  ggplot(aes(prop_monitored_lab, diff_passes)) +
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf, fill = colors[1], alpha = .12) + 
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0, fill = colors[2], alpha = .12) + 
  geom_text(data = data_rel, aes(label = label, color = color), hjust = 0, fontface = "bold", size = 6) +
  geom_hline(yintercept = 0, size = .6, linetype = "31", color = "grey75") +
  geom_boxplot(color = "grey45", width = .85, size = .5, 
               outlier.size = .8, outlier.alpha = .4, outlier.shape = 1) +
  stat_summary(geom = "point", shape = 18, size = 3, color = "#212121") +
  facet_wrap(~distribution, nrow = 2, dir = "v") + 
  scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  scale_y_continuous(expand = c(.03, .03), breaks = seq(-1, 3.5, by = .5)) +
  scale_color_identity() +
  labs(x = "Area covered by AUD",
       y = "Relative difference of predicted versus expected bat passes") +
  theme(axis.text.x = element_text(size = 11))
Show code
ggsave(here("plots", "figS3_passes_difference_boxplots_rel.pdf"), 
       width = 15, height = 8.2, device = cairo_pdf)
ggsave(here("plots", "figS3_passes_difference_boxplots_rel.png"), 
       width = 15, height = 8.2, dpi = 800)

Visualizations Fatalities

Show code
fits_scenarios <- 
  scenarios %>% 
  mutate(n_expected = ifelse(n_monitored == 0, 0, n_monitored / prop_monitored)) %>% 
  nest(data = -c(distribution, prop_monitored, prop_monitored_lab)) %>% 
  mutate(
    fit = map(data, ~ lm(n_fatalities ~ n_expected , data = .x)),
    tidied = map(fit, tidy),
    glance = map(fit, glance),
    augment = map(fit, augment),
    rsq = glance %>% map_dbl('r.squared'),
    slope = tidied %>% map_dbl(function(x) x$estimate[2]),
    se = tidied %>% map_dbl(function(x) x$std.error[2])
  ) %>% 
  dplyr::select(distribution, prop_monitored, prop_monitored_lab, rsq, slope, se) %>% 
  mutate(
    t = (.01 - slope) / se, ## t value = (expected slope - observed slope) / standard error
    p = 2*pt(abs(t), 1200-2, lower.tail = FALSE) ## p-value = two-tailed test of t value and d.f. (n - 2)
  )

Regression Fatalities vs Predicted Bat Passes — What we observe

Dedicated Scenarios

Low-Coverage Scenario
Show code
p_data_low <- 
  scenarios %>% 
  filter(prop_monitored == .04) %>% 
  mutate(n_expected = ifelse(n_monitored == 0, 0, n_monitored / prop_monitored))

p_data_low %>% 
  left_join(fits_scenarios) %>% 
  ggplot(aes(n_expected, n_fatalities)) +
  geom_abline(intercept = 0, slope = .01, color = "grey45", linetype = "31", size = .4) +
  geom_point(shape = 16, alpha = .08) +
  geom_quantile(quantiles = c(0.25, 0.75), color = "#D89684", alpha = .7, size = .6) +
  geom_quantile(quantiles = c(0.5), color = "#E60000", size = .9) +
  # geom_text(
  #   aes(label = paste("R² = ", sprintf("%1.2f", rsq)), 
  #       x = 0, y = 87), 
  #   stat = "unique", family = "Open Sans", color = "grey30",
  #   size = 3.7, hjust = 0
  # ) +
  facet_wrap(~ distribution, nrow = 2, scales = "free_x", dir = "v") +
  scale_x_continuous(breaks = 0:5*4000, limits = c(NA, max(p_data_low$n_expected)),
                     labels = scales::comma_format()) +
  scale_y_continuous(breaks = 0:4*20) +
  labs(title = "Low coverage scenarios",
       subtitle = "4% acoustic monitoring area",
       x = "Number of bat passes (predicted)", y = "Number of fatalities (simulated)") + 
  theme(
    panel.grid.major.y = element_blank(), 
    panel.spacing.x = unit(1.2, "lines"), 
    panel.spacing.y = unit(.6, "lines"), 
    legend.text = element_text(hjust = 1),
    strip.text = element_text(hjust = 0),
    plot.subtitle = element_text(margin = margin(t = 0, b = 20))
  )
Show code
ggsave(here("plots", "fig5_fatalities_correlation_low_coverage.pdf"), 
       width = 14, height = 7.5, device = cairo_pdf)

ggsave(here("plots", "fig5_fatalities_correlation_low_coverage.png"), 
       width = 14, height = 7.5, dpi = 800)

High-Coverage Scenario

Show code
p_data_high <- 
  scenarios %>% 
  filter(prop_monitored == .5) %>% 
  mutate(n_expected = ifelse(n_monitored == 0, 0, n_monitored / prop_monitored))

p_data_high %>% 
  left_join(fits_scenarios) %>% 
  ggplot(aes(n_expected, n_fatalities)) +
  geom_abline(intercept = 0, slope = .01, color = "grey45", linetype = "31", size = .4) +
  geom_point(shape = 16, alpha = .08) +
  geom_quantile(quantiles = c(0.25, 0.75), color = "#D89684", alpha = .7, size = .6) +
  geom_quantile(quantiles = c(0.5), color = "#E60000", size = .9) +
  # geom_text(
  #   aes(label = paste("R² = ", sprintf("%1.2f", rsq)), 
  #       x = 0, y = 87), 
  #   stat = "unique", family = "Open Sans", color = "grey30",
  #   size = 3.7, hjust = 0
  # ) +
  facet_wrap(~ distribution, nrow = 2, scales = "free_x", dir = "v") +
  scale_x_continuous(breaks = 0:5*2000, limits = c(NA, max(p_data_high$n_expected)),
                     labels = scales::comma_format()) +
  scale_y_continuous(breaks = 0:4*20) +
  labs(title = "High coverage scenarios",
       subtitle = "50% acoustic monitoring",
       x = "Number of bat passes (predicted)", y = "Number of fatalities (simulated)") + 
  theme(
    panel.grid.major.y = element_blank(), 
    panel.spacing.x = unit(1.2, "lines"), 
    panel.spacing.y = unit(.6, "lines"), 
    legend.text = element_text(hjust = 1),
    strip.text = element_text(hjust = 0),
    plot.subtitle = element_text(margin = margin(t = 0, b = 20))
  )
Show code
ggsave(here("plots", "fig6_fatalities_correlation_high_coverage.pdf"), 
       width = 14, height = 7.5, device = cairo_pdf)

ggsave(here("plots", "fig6_fatalities_correlation_high_coverage.png"), 
       width = 14, height = 7.5, dpi = 800)

All Scenarios with 400 Bat Passes

Show code
p_data_all <- scenarios %>% 
  mutate(n_expected = ifelse(n_monitored == 0, 0, n_monitored / prop_monitored)) 

p_data_all %>% 
  ggplot(aes(n_expected, n_fatalities)) +
  geom_abline(intercept = 0, slope = .01, color = "grey45", linetype = "31", size = .4) +
  geom_point(shape = 16, alpha = .08) +
  geom_quantile(quantiles = c(0.25, 0.75), color = "#D89684", alpha = .7, size = .6) +
  geom_quantile(quantiles = c(0.5), color = "#E60000", size = .9) +
  # geom_text(
  #   aes(label = paste("R² = ", sprintf("%1.2f", rsq)), 
  #       x = 0, y = 87), 
  #   stat = "unique", family = "Open Sans", color = "grey30",
  #   size = 3.7, hjust = 0
  # ) +
  facet_grid(prop_monitored_lab ~ distribution, scales = "free_x") +
  scale_x_continuous(breaks = 0:5*4000, limits = c(NA, max(p_data_all$n_expected)),
                     labels = scales::comma_format()) +
  scale_y_continuous(breaks = 0:4*20) +
  labs(x = "Number of bat passes (predicted)", y = "Number of fatalities (simulated)") + 
  theme(
    panel.grid.major.y = element_blank(), 
    panel.spacing = unit(.8, "lines"), 
    legend.text = element_text(hjust = 1),
    strip.text.x = element_text(hjust = 0),
    strip.text.y = element_text(angle = 0, hjust = 1, margin = margin(3, 9, 3, 9))
  ) -> d

ggsave(here("plots", "figS4_fatalities_correlation_all.pdf"), 
       width = 24, height = 15, device = cairo_pdf)

ggsave(here("plots", "figS4_fatalities_correlation_all.png"), 
       width = 24, height = 15, dpi = 800)

p-Values Slopes

Show code
fits_scenarios %>% 
  ggplot(aes(prop_monitored_lab, fct_rev(distribution))) +
  #geom_tile(aes(fill = base), color = "white", size = .7) +
  geom_tile(aes(fill = p), color = "white", size = .8) +
  geom_text(
    aes(label = sprintf("%1.2f", p), color = p < .05),
    family = "Open Sans", size = 4.5
  ) +
  coord_cartesian(expand = FALSE, clip = "off") +
  scale_color_manual(values = c("#fefefe", "#212121"), guide = "none") +
  scale_x_discrete(position = "top") +
  # scico::scale_fill_scico(
  #   palette = "batlowW", direction = -1, name = "P-value (t-test, df = 1.198)",
  #   begin = .12, breaks = seq(0, 1, by = .1), midpoint = .05, limits = c(0, 1)
  # ) +
  scico::scale_fill_scico(
    palette = "batlowW", direction = -1, end = .85, 
    name = "P-value (t-test, df = 1.198)", limits = c(0, 1),
    breaks = seq(0, 1, by = .1), labels = function(x) sprintf("%1.1f", x)
  ) +
  scale_alpha_manual(values = c(0, 1), guide = "none") +
  scale_alpha_manual(values = c(0, 1), guide = "none") +
  guides(fill = guide_colorsteps(title.position = "top", show.limits = TRUE)) +
  labs(x = "Proportion covered by AUD", y = NULL) +
  theme(panel.spacing = unit(.9, "lines"), 
        panel.background = element_rect(size = .7, color = "white", fill = "transparent"),
        panel.border = element_rect(color = "transparent", fill = "transparent"),
        axis.title.x.top = element_text(face = "bold", margin = margin(t = 8, b = 12)),
        axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14, face = "bold", hjust = 0),
        axis.ticks = element_line(color = "transparent"),
        axis.ticks.length = unit(.3, "lines"),
        legend.position = "top",
        legend.justification = "left",
        legend.key.width = unit(6, "lines"), legend.key.height = unit(.6, "lines"))
Show code
ggsave(here("plots", "tab1_slopes_ttest_heatmap.pdf"), 
       width = 12, height = 6.7, device = cairo_pdf)
ggsave(here("plots", "tab1_slopes_ttest_heatmap.png"), 
       width = 12, height = 6.7, dpi = 800)

Additional Visualizations

Heatmap Deviation Bat Passes Recorded vs. Expected

Boxplots Proportion of Passes Monitored

Show code
box <- 
  ggplot(scenarios, aes(factor(n), prop_n_monitored)) +
  geom_hline(aes(yintercept = prop_monitored), color = "grey75", size = .7)  + 
  geom_boxplot(aes(color = diff), size = .4, outlier.shape = 1, outlier.size = .3) + 
  scale_y_continuous(labels = scales::percent_format(), 
                     sec.axis = dup_axis(name = "Area covered by AUD", breaks = NULL, labels = NULL)) +
  scale_color_gradient2(low = "firebrick", mid = "grey40", high = "firebrick", 
                        limits = c(-.5, .5), guide = "none") +
  labs(x = "Number of bat passes", y = "Proportion of passes monitored") +
  theme(panel.spacing.y = unit(.65, "lines"), axis.text = element_text(size = 9),
        axis.title.x = element_text(size = 18, face = "bold"), 
        axis.title.y = element_text(size = 18, face = "bold"), 
        strip.text.x = element_text(size = 18),
        strip.text.y = element_text(size = 16),
        axis.title.y.right = element_text(size = 21, color = "grey60", face = "bold", margin = margin(l = 12)))

## fixed y scale
box + facet_grid(prop_monitored_lab ~ distribution)
Show code
ggsave(here("plots", "bat_passes", "passes_recorded_boxplots_fixed.pdf"), 
       width = 23, height = 13, device = cairo_pdf)
ggsave(here("plots", "bat_passes", "passes_recorded_boxplots_fixed.png"), 
       width = 23, height = 13, dpi = 800)

## free y scale
box + facet_grid(prop_monitored_lab ~ distribution, scale = "free_y") 
Show code
ggsave(here("plots", "bat_passes", "passes_recorded_boxplots_free.pdf"), 
       width = 23, height = 13, device = cairo_pdf)
ggsave(here("plots", "bat_passes", "passes_recorded_boxplots_free.png"), 
       width = 23, height = 13, dpi = 800)

Boxplots Absolute Difference Predicted vs. Expected

for dedicated scenarios

Show code
data_abs_scen <- 
  scenarios %>% 
  filter((n == 400 & prop_monitored == .04) | (n == 100 & prop_monitored == .5)) %>% 
  mutate(
    diff_passes = (n_monitored / prop_monitored) - n,
    scenario = if_else(prop_monitored == .04, "Low", "High")
  ) %>% 
  bind_rows(
    tibble(
      distribution = factor("a) uniform", levels = levels(scenarios$distribution)),
      label = c("Overestimation", "Underestimation"),
      scenario = "High",
      diff_passes = c(625, -250),
      color = colors
    )
  )
  
data_abs_scen %>%
  filter(is.na(label)) %>% 
  ggplot(aes(scenario, diff_passes)) +
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf, fill = colors[1], alpha = .12) + 
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0, fill = colors[2], alpha = .12) + 
  geom_text(data = data_abs_scen, aes(label = label, color = color), hjust = .15, fontface = "bold", size = 5) +
  geom_hline(yintercept = 0, size = .6, linetype = "31", color = "grey75") +
  geom_boxplot(color = "grey45", size = .6, width = .5, 
               position = position_nudge(x = .05), 
               outlier.size = .9, outlier.shape = 1) +
  geom_point(color = "grey90", shape = "-", size = 10.1, 
             position = position_nudge(x = -.31)) +
  geom_point(color = "grey45", shape = "-", size = 10.1, alpha = .05, 
             position = position_nudge(x = -.31)) +
  stat_summary(geom = "point", shape = 18, size = 4.5, color = "#212121", 
               position = position_nudge(x = .05)) +
  facet_wrap(~distribution, nrow = 2, dir = "v") + 
  scale_y_continuous(expand = c(.04, .04), breaks = seq(-800, 1600, by = 400)) +
  scale_color_identity() +
  labs(x = "Coverage scenario",
       y = "Absolute difference of predicted versus expected bat passes") +
  theme(axis.text.x = element_text(size = 13))
Show code
ggsave(here("plots", "bat_passes", "passes_difference_boxplots_abs_scenarios.pdf"), 
       width = 12, height = 7.7, device = cairo_pdf)
ggsave(here("plots", "bat_passes", "passes_difference_boxplots_abs_scenarios.png"), 
       width = 12, height = 7.7, dpi = 800)

All Scenarios with 400 Bat Passes

Show code
## absolute
data_abs <- 
  scenarios %>%
  filter(n == 400) %>% 
  filter(!prop_monitored == .04) %>% ## remove dedicated scenario 
  mutate(diff_passes = (n_monitored / prop_monitored) - n) %>% 
  bind_rows(
    tibble(
      distribution = factor("a) uniform", levels = levels(scenarios$distribution)),
      label = c("Overestimation", "Underestimation"),
      prop_monitored_lab = factor("10%", levels = levels(scenarios$prop_monitored_lab)), 
      diff_passes = c(625, -250),
      color = colors
    )
  )

data_abs %>% 
  filter(is.na(label)) %>% 
  ggplot(aes(prop_monitored_lab, diff_passes)) +
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf, fill = colors[1], alpha = .12) + 
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0, fill = colors[2], alpha = .12) + 
  geom_text(data = data_abs, aes(label = label, color = color), hjust = 0, fontface = "bold", size = 6) +
  geom_hline(yintercept = 0, size = .6, linetype = "31", color = "grey75") +
  geom_boxplot(color = "grey45", width = .85, size = .5, 
               outlier.size = .8, outlier.alpha = .4, outlier.shape = 1) +
  stat_summary(geom = "point", shape = 18, size = 3, color = "#212121") +
  facet_wrap(~distribution, nrow = 2, dir = "v") + 
  scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  scale_y_continuous(expand = c(.03, .03), breaks = seq(-800, 1600, by = 400)) +
  scale_color_identity() +
  labs(x = "Area covered by AUD",
       y = "Absolute difference of predicted versus expected bat passes") +
  theme(axis.text.x = element_text(size = 11))
Show code
ggsave(here("plots", "bat_passes", "passes_difference_boxplots_abs.pdf"), 
       width = 15, height = 8.2, device = cairo_pdf)
ggsave(here("plots", "bat_passes", "passes_difference_boxplots_abs.png"), 
       width = 15, height = 8.2, dpi = 800)

Regression Fatalities vs Predicted Bat Passes — What really happens

Show code
## with variation
scenarios %>% 
  filter(prop_monitored == .05) %>% 
  ggplot(aes(n, n_fatalities)) +
  #geom_abline(slope = .01, intercept = 0, size = .6, color = "grey75") +
  geom_jitter(alpha = .02, width = 20, height = 0) +
  geom_quantile(quantiles = c(.25, .75), color = "dodgerblue3", size = .6) +
  geom_quantile(quantiles = c(.5), color = "red", size = .9) +
  facet_grid(~ distribution) +
  scale_x_continuous(expand = c(.1, .1)) +
  labs(x = "Number of bat passes (reality)",
       y = "Number of fatalities") +
  theme(panel.grid.major.y = element_blank(), axis.text.x = element_text(size = 10))
Show code
ggsave(here("plots", "fatalities", "fatalities_correlation_reality_var.pdf"), 
       width = 13.5, height = 3.2, device = cairo_pdf)
ggsave(here("plots", "fatalities", "fatalities_correlation_reality_var.png"), 
       width = 13.5, height = 3.2, dpi = 800)


## without variation
# scenarios %>% 
#   filter(prop_monitored == .05) %>% 
#   mutate(n_fatalities = round(n / 100)) %>% 
#   ggplot(aes(n, n_fatalities)) +
#   #geom_abline(slope = .01, intercept = 0, size = .6, color = "grey75") +
#   geom_point(alpha = .02) +
#   geom_quantile(quantiles = c(0.25, 0.75), color = "dodgerblue3", size = .6) +
#   geom_quantile(quantiles = c(0.5), color = "red", size = .9) +
#   facet_grid(~ distribution) +
#   scale_x_continuous(expand = c(.1, .1)) +
#   labs(x = "Number of bat passes (reality)",
#        y = "Number of fatalities") + #,
#        #title = "What really happens",
#        #subtitle = "Number of fatalities per number of bat passes (1% correlation without variation)") +
#   theme(panel.grid.major.y = element_blank(), axis.text.x = element_text(size = 10))
# 
# ggsave(here("plots", "fig_Sx_fatalities_correlation_reality_fix.pdf"), 
#        width = 13.5, height = 3.2, device = cairo_pdf)
# ggsave(here("plots", "fig_Sx_fatalities_correlation_reality_fix.png"), 
#        width = 13.5, height = 3.2, dpi = 800)

Session Info
Show code
[1] "2022-03-31 17:02:22 CEST"
Show code
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    
system code page: 65001

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
 [1] patchwork_1.1.1 broom_0.7.11    here_1.0.1      forcats_0.5.1  
 [5] stringr_1.4.0   dplyr_1.0.7     purrr_0.3.4     readr_2.0.2    
 [9] tidyr_1.1.4     tibble_3.1.6    ggplot2_3.3.5   tidyverse_1.3.1

loaded via a namespace (and not attached):
  [1] colorspace_2.0-2     ellipsis_0.3.2       class_7.3-19        
  [4] rprojroot_2.0.2      fs_1.5.0             rstudioapi_0.13     
  [7] listenv_0.8.0        farver_2.1.0         MatrixModels_0.5-0  
 [10] scico_1.3.0          prodlim_2019.11.13   fansi_0.5.0         
 [13] lubridate_1.8.0      xml2_1.3.2           codetools_0.2-18    
 [16] splines_4.1.2        downlit_0.4.0        cachem_1.0.6        
 [19] knitr_1.36           jsonlite_1.7.2       pROC_1.18.0         
 [22] caret_6.0-90         dbplyr_2.1.1         compiler_4.1.2      
 [25] httr_1.4.2           backports_1.2.1      assertthat_0.2.1    
 [28] Matrix_1.3-4         fastmap_1.1.0        cli_3.1.0           
 [31] htmltools_0.5.2      quantreg_5.86        tools_4.1.2         
 [34] gtable_0.3.0         glue_1.4.2           reshape2_1.4.4      
 [37] Rcpp_1.0.7           cellranger_1.1.0     jquerylib_0.1.4     
 [40] vctrs_0.3.8          nlme_3.1-153         conquer_1.2.1       
 [43] iterators_1.0.13     timeDate_3043.102    xfun_0.27           
 [46] gower_0.2.2          globals_0.14.0       rvest_1.0.2         
 [49] lifecycle_1.0.1      future_1.22.1        MASS_7.3-54         
 [52] scales_1.1.1         ipred_0.9-12         ragg_1.1.3          
 [55] hms_1.1.1            parallel_4.1.2       SparseM_1.81        
 [58] yaml_2.2.1           memoise_2.0.1        sass_0.4.0          
 [61] distill_1.3          rpart_4.1-15         stringi_1.7.5       
 [64] highr_0.9            foreach_1.5.1        lava_1.6.10         
 [67] rlang_0.4.12         pkgconfig_2.0.3      systemfonts_1.0.3   
 [70] matrixStats_0.61.0   evaluate_0.14        lattice_0.20-45     
 [73] labeling_0.4.2       recipes_0.1.17       tidyselect_1.1.1    
 [76] parallelly_1.29.0    plyr_1.8.6           magrittr_2.0.1      
 [79] R6_2.5.1             generics_0.1.1       DBI_1.1.2           
 [82] pillar_1.6.4         haven_2.4.3          withr_2.4.3         
 [85] survival_3.2-13      nnet_7.3-16          future.apply_1.8.1  
 [88] modelr_0.1.8         crayon_1.4.2         utf8_1.2.2          
 [91] tzdb_0.1.2           rmarkdown_2.11       grid_4.1.2          
 [94] readxl_1.3.1         data.table_1.14.2    ModelMetrics_1.2.2.2
 [97] reprex_2.0.1         digest_0.6.29        textshaping_0.3.6   
[100] stats4_4.1.2         munsell_0.5.0        bslib_0.3.1